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We present tests and results of a new axisymmetric, fully general relativistic code 
capable of solving the coupled Einstein-matter system for a perfect fluid matter 
field. Our implementation is based on the Bondi metric, by which the spacetime 
is foliated with a family of outgoing light cones. We use high-resolution shock- 
capturing schemes to solve the fluid equations. The code can accurately maintain 
long-term stability of a spherically symmetric, relativistic, polytropic equilibrium 
model of a neutron star. In axisymmetry, we demonstrate global energy conserva- 
tion of a perturbed neutron star in a compactified spacetime, for which the total 
energy radiated away by gravitational waves corresponds to a significant fraction 
of the Bondi mass. 



1 Characteristic numerical relativity and hydrodynamics 

1.1 Introduction 

We describe a new axisymmetric, fully general relativistic numerical code evolving 
tiie Einstein equations along with perfect fluid matter and apply it to studies of 
neutron stars. Our numerical implementation of theJield equations of general rel- 
ativity is based on the light cone formalism of BondiH and Tamburino-WinicouiO. 
The light cone approach has a number of advantages compared to spacelike folia- 
tions, i) It is physically motivated; the light cones offer a simple and unambiguous 
physical gauge on which to base the numerical spacetime grid, ii) It is uncon- 
strained; the evolved variables capture rather directly the true degrees of freedom 
of the gravitational field, iii) It is very efScient; there are but two partial differ- 
ential equations to solve, along with a set of radial integrations, iv) It allows for 
well defined compactification of the domain, which leads to perfect outer boundary 
conditions, v) Finally, and may be most importantly, the above theoretical advan- 
tages have been shown in a series of worJts to translate to remarkably numerically 
robust aruLstable codes, see for example 13. For a recent review of the approach see 
Winicourtj. 

The broad target of this project is the detailed investigation of neutron star 
dynamics using numerical relativity. A prerequisite for such studies is the devel- 
opment of very accurate and long-term stable general relativistic codes. It would 
seem that the feature list of the characteristic approach makes it ideal for such 
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studies. One serious bottleneck of the approach immediately recognized by the 
trained relativist is the breakdown of a lightlike coordinate system in the emer- 
gence of light caustics. Interestingly though, due to its quasispherical nature, no 
matter how agitated a neutron star, it is unlikely to focus the light cones emanat- 
ing from its interior. A limitation of more technical nature is the existence of a 
somewhat restrictive Courant-Friedrichs-Lewy-condition for explicit algorithms in 
multidimensions. Implicit evolution schemes would help in this direction. 

The incorporation of fluid matter fields in the |Characteristic formulation of the 
Einstein equation was considered as early as 19833 and numerical results were re- 
ported in referencetJ. A different line of attack was initiated recentlytll which brings 
into the considerations the modern machinery from Computational Fluid Mechan- 
ics. In this approach, the evolution equations for the lamtter fields are solved us- 
ing relativistic high-resolution shock-capturing schemesoO (HRSC schemes) based 
upon (exact or approximate) Riemann solvers. At this early stage of our investiga- 
tions, the ability of capturing shock waves in the hydrodynamics is not essential. 
Nevertheless it will be so when studying core-collapse scenarios. The implementa- 
tion of HRSC schemes in a general characteristic code without imposing symmetry 
conditions is the current subject of a collaboration (GRACE, General relativistic 
astrophysics via characteristic evolution). 

In this work we restrict ourselves to axisymmetric spacetimes. The code has 
been developed building on previous work by Gomez, Papadopoulos and WinicourQ 
(GPW), which constructed an axisymmetric characteristic vacuum code, to which 
we have added a perfect fluid matter field. ApplicationH in spherical, dynamic black 
hole spacetimes have been presented in referenceaiil'E£l. Studies of_the spherically 
symmetric collapse of supermassive stars are discussed in referencdj. 



1.2 Fundamental equations 

The main equations we have to solve for the problem we are aiming at are the 
Einstein equation 

Gab - ^^Tab = (1) 

and the conservation equation of the stress-energy tensor Tab, 

V^Tab = 0. (2) 

The stress-energy tensor is chosen as that of a perfect fluid. 

Tab = phUaUb + P9ab, (3) 

where p denotes the mass density, h the relativistic specific enthalpy and Ua the 
four-velocity of the fluid. The pressure p is subject to an equation of state. We 
choose a perfect fluid equation of state 

p={V-l)pe, (4) 

where e denotes the specific internal energy of the fluid, related to the specific 
enthalpy asft, = l + e+ ^. In addition we assume conservation of baryonic matter, 

V"(pu,) = 0. (5) 
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We solve the above equations using the standard coordinate system of characteristic 
numerical relativity (see Figure |l[). 




Figure 1. Foliation of the spaeetime by null hypersurfaces (light cones): Starting from the geodesic 
of a freely falling particle, which defines the origin of the coordinate system, surface forming light 
rays are emitted. We define a new time coordinate u to be constant along each outgoing light 
cone. The Killing coordinate {j> is suppressed in the diagram. The origin of the coordinate system 
will lie (in the axisymmetric case) on the symmetry axis of the star. If the system has additional 
reflection symmetry with respect to the equator, the origin is actually at the center of the star. 



More precisely, we use the Bondi metric 

ds^ = -(-6^/3 _ U\^e^^)du^ - 2e^^du dr - 2Ur^e^^du d0 
r 

+r^{e^^de^ + e-^"'sin^0 dcj)^) (6) 

with null coordinate u, radial coordinate r, azimuthal coordinate 9 and the spher- 
ical coordinate (p, which is a Killing coordinate. The four metric fields /3, 7, U 
and V are determined by solving the Einstein equation (Q). With a choice of a 
coordinate chart (cc^jX*), the hydrodynamic equations (j^) and are written as 
an explicitJuiperbolic system, well-suited for numerical applications, as detailed in 
referenceal2l't3. 



1.3 Numerical implementation 

The numerical implementation of the Einstein equation (|l|) closely follows that of 
GPW. Using the Bondi metric (^, equation (|l|) splits into hypersurface equations on 
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each light cone (for the fields (3, U and V) and one evolution equation (for 7), a wave 
equation. These equations are solved by the same marching algorithms described 
in GPW, now with additional source terms arising from the matter fields. As 
we are using spherical coordinates, we have to take special care with the numerical 
treatment of the coordinate singularities at the origin and the polar axis. To this aim 
we extend the treatment of the origin presented in GPW to nonvacuum spacetimes 
and regularize the poles by introducing a new azimuthal coordinate y = —cos9 
and by redefining variables, for example 7 = 7 sin^d. As we want to keep the 
freedom of working with numerical grids which only cover the neutron star without 
its vacuum exterior, we generalize the radial coordinate used in GPW. Starting 
from an equidistant radial coordinate x G [0, 1], we allow for a general coordinate 
transformation of the form r ~ r{x), which enables us to use compactified or non- 
compactified grids. 

As already mentioned, we use HRSC schemes to solve the fluid equations (||) 
and (H). They constitute a hyperbolic system of balance laws. We use a so-called 
method of lines written in conservation form to solve these equations and to update 
the initial data forward in time. The principal part is solved by a Giidunov-like 
method, with the same approximate Riemann solver as in referencetil (see this 
reference for further details). 

2 Tests and Applications 

In this section, we apply our coupled code to models of neutron stars. For com- 
pleteness we mention that the axisymmetric vacuum part of the code has been 
successfully tested in additional situations. These comprise the comparison to an 
exact solution (SIMPLE), the evolution of weak ingoing gravitational wavea,and an 
energy conservation test for vacuum data, following the tests described in 

2.1 Long-term stability of spherical neutron stars 

As a simplified model for a self-gravitating neutron star we consider the spherically 
symmetric solution of the general relativistic hydrostatic equation with a polytropic 
equation of state, p — K/f^- This is the Tolman-Oppenheimer-Volkoff solution in its 
light cone representation^. In our analysis, we have wopked with two representative 
models, following the choice of parameters made also in a. We restrict the discussion 
to the more relativistic one of these models. In dimensionless units {c — G ~ 
Mq — 1), the equilibrium properties of this star are described by: a polytropic 
index F = 2, a polytropic constant K = 100 and a central density pc = 1.28 lO^'^. 
The equilibrium model has a total mass M = 1.4. The time light needs to travel 
across this neutron star corresponds to about 17u in our time unit. We use this 
model for convergence tests as well as to check long-term stability of the numerical 
algorithms. 

As the Tolman-Oppenheimer-Volkoff solution is static, convergence can be easily 
checked by subtracting the evolved solution from the initial (true) solution. We 
find second order convergence in all variables. Figure ^ shows a suitable norm of 
all variables as a function of the radial coordinate, which measures the deviation 
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from the initial solution for different grid resolutions. 
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Figure 2. Convergence test for the Tolmann-Oppenheimer-VolkofF solution. Plotted is a suitable 
norm of all variables, which measures the deviation from equilibrium, for the indicated number of 
radial grid points (gp), after about two light-crossing times. To good approximation the solution 
is 2nd order convergent. 

Figure || displays the L2-norm of the radial velocity ||w''||2 of the star as a 
function of time. Due to discretization errors radial fluid motion is generated, the 
radial velocity slightly deviates from zero and the star oscillates in its radial modes 
of pulsation. These radial oscillations do not increase with time, which reflects the 
long-term stability of our numerical implementation. 

Figure ^ shows the density profile of the neutron star for a very long integration 
time u=10000 for a lower resolution. Even though this corresponds to a very long- 
term, fully general relativistic hydrodynamic evolution, the density profile almost 
does not change, the equilibrium of the star being maintained to very good precision. 
The upper line corresponds to the initial model, the lowest line to the density profile 
at the time u=10000, which corresponds to roughly 590 light-crossing times. 

2.2 Non-spherical matter motions 

Next we consider the spherical model of the neutron star described in the previous 
section and add a non-spherical density component. This will lead to non-spherical 
evolution and hence allows for testing the non-spherical terms of the evolution sys- 
tem. More precisely, we add a Gaussian profile, which is centered at half radius at 
the polar axis and whose maximal value corresponds to 10% of the density of the 
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Figure 3. Radial velocity of the star, averaged over the radial coordinate, as a function of time. 
The evolution corresponds to roughly 59 light-crossing times. 
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Figure 4. Density profile of the neutron star at the initial and subsequent times. The final time 
corresponds to u = 10000 (roughly 590 light-crossing times). The equilibrium model is maintained 
to a very good precision after such very long integration times. 
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equilibrium model at that point. Our initial model thus shows strong deviations 
from spherical symmetry. The non-spherical fluid motion will trigger gravitational 
waves. Figure || shows a series of snapshots of this evolution. We plot here the 
evolution of the overdensity, that is, the difference between the non-spherical so- 
lution and the equilibrium one. The evolution induces a large mass transfer. Due 
to this mass transfer, we would expect a remarkable gravitational wave signal, i.e. 
the quantity 7 picks up a non-zero contribution. We found that our code is also 
second order convergent even for such non-spherical data. As we are in a regime 
where the exact solution is not known, we use Cauchy convergence for this test. 




Figure 5. Evolution of a non-spherical density distribution on top of the equilibrium model of the 
neutron star. Plotted are the profile of the initial overdensity and the profiles of the overdensity 
at subsequent times as a function of the coordinates r and y. The evolution shows a large, non- 
spherical mass transfer, which allows us to check the second-order accurate implementation of all 
terms in our code. 
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2.3 Generation of gravitational waves 

One of the major benefits of using lightlike foliations is the possibility of correctly 
extracting gravitational radiation after a compactification of spacetime. For an 
asymptotically flat spacetime, the metric fields can be Taylor expanded around fu- 
ture null infinity 2""*" (which corresponds to the limit r — > cxd using our null foliation) 
in inverse powers of the radial coordinate r. Using the coefficients of this power 
series, the radiated energy emitted by gravitational waves can be determined.Q In 
the same way, one could determine the Bondi mass, the total energy. However, as 
one would have to pick up nonleading terms in the power series expansion for the 
metric variable V , as the leading terms diverge on X"*", the numerical extraction is 
easily spoiled by numerical errors. We therefore follow the work of Gomez et aid 
and globally introduce new metric variables, from which the Bondi mass can be 
determined using the leading term of their power series expansion. 

In this section, we focus on a global energy conservation test. We start with a 
strongly perturbed neutron star and evolve this configuration for a very short time, 
which is sufficient for our purpose of testing convergence. We use a strong ingoing 
gravitational wave to perturb the equilibrium star, 

^ = 0.05e-2('-4)%-4^'. (7) 

Such a large amplitude is not realistic, but we choose it to test our numerical 
implementation in the nonlinear regime. Such initial data involve large velocities in 
the fluid of the star and produce strong outgoing gravitational waves. By 'strong' we 
understand that the energy which is radiated away by gravitational waves is larger 
than the numerical errors for the calculation of the Bondi mass for the chosen 
resolution and integration time. Let M be the Bondi mass and P the total energy 
radiated away by gravitational waves. Thus, the convergence of the quantity 

ec := M\u=a - M\u=u,>o - ^'l[o,u.] (8) 

to zero represents a very severe global test for our implementation (see Figure ^ . 

The obtained first-order convergence rate can be explained by the use of a total 
variation diminishing HRSC scheme, which, although it is second-order accurate in 
smooth, monotonous parts of the fiow, reduces to first-order at local extrema, which 
are present in the interior of the numerical domain in this test. Tests including 
propagation and scattering off the origin of pure vacuum gravitational fields yielded 
the expected second-order convergence. 

There are_two additional consistency conditions, which relate the metric quan- 
tities on X+ .Ej These are 

S ^U,e + U cote (9) 
7. = -r-^^^^neie^^^^U (10) 



where S is globally defined as 



V -r 

s^—^- (11) 



Here again, we find first-order convergence. 
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Figure 6. Global energy conservation test for a neutron star and a strong gravitational wave. 
Plotted are the deviation from energy conservation (triangles) and the total energy emitted by 
gravitational waves (squares). The final integration time is ti* = 0.002, the number of angular 
grid points is 0.2 gp + 1, where gp denotes the number of radial grid points. 

3 Conclusion 

We have developed an axisymmetric, fully general relativistic characteristic code 
with a perfect fluid matter fleld. We have checked convergence in spacetimes con- 
taining a neutron star. For a spherical, relativistic equilibrium model of a poly- 
tropic neutron star, we have shown long-term stability for hundreds of light-crossing 
times. In axisymmetry, we analyzed convergence in strongly perturbed neutron 
stars. Compactifying the spacetime, we could show global energy conservation, 
in the sense that the amount of energy lost by the system is radiated away by 
gravitational waves. 
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